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Temperature measurement under the convection and 
segregation in the vibrated bed of powder: A numerical study 
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Abstract In numerically simulated vibrated beds of powder, 
we measure temperature under convection by the generalized 
Einstein's relation. The spatial temperature distribution turns 
out to be quite uniform except for the boundary layers. In ad- 
dition to this, temperature remains uniform even if segregation 
occurs. This suggests the possibility that there exists some ther- 
mal equilibrium state even in a vibrated bed of powder. This 
finding may lead to a unified view of the dynamic steady state 
of granular matter. 
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Introduction 

Granular matter is characterized by a set of macroscopic parti- 
cles, e.g., sand, snow fall, sugar, salt and rice. Recently many 
physicists became interested in the dynamics of granular mat- 
ter [1 1, because it shows many striking features. These include 
convection, segregation, surface waves, and flows. They are 
hard to understand by the word of physics, although it is rel- 
atively easy to reproduce them numerically. 

For example, the concept of granular temperature [ 1 1 is of- 
ten used to describe the dynamical state of powder. However, 
in contrast to rapid granular flow [2|, the granular temperature 
turns out not to be a good measure for dense granular matter. 
For example, granular temperature in vibrated beds of powder 
is anisotropic |3], which is not a favorable feature. 

On the other hand, Wildman and Huntley 1 4 1 have found 
useful a temperature that is given by the mean squared dis- 
placement curves for short times. Although it sounds promis- 
ing, their definition is valid only when convective flow is neg- 
ligible. 

Recently Makse and Kurchan 1 5 1 have numerically found 
that the temperature defined by Einstein's relation is well de- 
fined and is independent of the size of the particles in gran- 
ular matter. In contrast to that of Wildman and Huntley |4|, 
Makse and Kurchan's definition is valid even if there is a non- 
negligible flow. Also D'Anna et al |6| have found that their 
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variant using AC frequency is also well defined. Thus we may 
expect that transport phenomena of granular particles can give 
us a more suitable definition of temperature than the conven- 
tional granular temperature. 

In spite of these successful applications of Einstein's rela- 
tion to define the temperature in the granular matter, its applica- 
tion is very limited. This is because the conventional Einstein's 
relation can be defined only in a spatially uniform system. In 
order to apply Einstein's relation, we need to understand the 
long time behavior of tagged particles. If the system is not 
uniform at all, the tagged particle will travel around various 
places having various temperatures. Thus what we can mea- 
sure is a spatially averaged temperature, which cannot reflect 
spatial structure in the system. On the other hand, if the system 
is uniform, temperature is uniform, too. This means that we 
can never know if the thermal equilibrium really stands or not 
from the measurement of temperature by Einstein's relation. 
However, as demonstrated by Wildman and Huntley |4), it is 
not always necessary to observe long time behavior in order to 
derive temperature from the transport phenomena. 

In this paper, we have developed a heuristic method to mea- 
sure temperature by Einstein's relation in the strongly non- 
uniform system. In [|2] following Makse and Kurchan |5|, we 
introduce the temperature using Einstein's relation. In the 
heuristic approach to define local temperature is presented, and 
in <j4]this definition is further generalized so as to measure tem- 
perature in the system under the flow. The numerical results 
in vibrated beds are shown in f|5] and how segregation affects 
the temperature is investigated in Sj6] In fJT] we have discussed 
about the validity of our results; fj8]contains summary and con- 
clusion. 



Einstein's relation and temperature 

Suppose there is a particle which obeys a random walk. So as to 
be simple, we restrict its motion to the one dimensional space. 
If x{t) is the coordinate of the particle at time t, then the diffu- 
sion constant D is defined as 

({x(t)-x(0)] 2 )=2Dt, 

where (• • •) means the ensemble average. On the other hand, 
the mobility /i can be defined as 

(x(t) - x(0)) = tft, 

where / is the applied force. When combining the above two 
equations with the Einstein's relation 

D 



T = 



2 



where T is temperature, we get 

f {[x(t)-x(0)f) 
2 (x(t) - x{0)) ' 

Recently, Makse and Kurchan |5| used this equation as the 
definition of temperature in granular matter. They have per- 
formed three dimensional distinct element method (DEM) for 
the periodic sheared system consisting of a binary set of both 
large and small particles. They measured D and /i along the di- 
rection perpendicular to the shear and found that the measured 
temperature is the same for both large and small particles. 

More recently, D' Anna et al | 6 1 have experimentally found 
that the fluctuation-dissipation ratio 

S{uS)oj 

can be a well defined temperature for the vibrated bed, where 
S(uj) is the noise power spectrum density of the angular fre- 
quency oj and x" is the imaginary part of the complex sus- 
ceptibility. Since this can be regarded as the generalization of 
the temperature introduced by Makse and Kurchan to the fre- 
quency dependent one, these kinds of definition of temperature 
seem to be valid for granular matter. 

However, there is a difficulty to apply the above definitions 
to general cases. This is because it is impossible to measure 
the spatial dependence of temperature using these definitions. 
Makse and Kurchan's definition is valid only for the long time 
limit, and D'Anna et a/'s procedure is used to measure tem- 
perature only when we can regard a whole vibrated bed as a 
thermal bath. Thus, it is rather difficult to relate temperature 
to the statistical mechanics, because in the statistical mechan- 
ics local temperatures must be defined. The definitions above 
cannot provide such information. 



Thus 

D = P out (i)£ 2 /2 = f P out (i)£ 2 
1 V J(i)/f 2 J(i) 

is the local temperature at cell i. Of course we do not insist 
that this is rigorous, but only regard it as a heuristic argument. 
Also, in order to remove the dependency on / we employ the 
following definition 

P eff = lim Ti 
of the effective temperature. 
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The definition of effective temperature in the system 
with flow 

As an example to test the above heuristic definition of the ef- 
fective temperature, we employ a numerical simulation of a 
vibrated bed of powder. In the vibrated bed of powder con- 
sidered here, there are convective motions of particles. Thus 
(x(t) — x(0)) is not zero even if / is equal to 0. Furthermore, 
the vibrated bed we deal with is a two dimensional system. 
Therefore we need some additional modification of the previ- 
ous definition of the effective temperature Tf . 

First we subdivide the system into N\ x N2 cells, where 
Ni£ is the horizontal size of system and N^l is the vertical size 
of the system. Then we denote the cell as (i, j ) if the cell is the 
ith in horizontal direction and the jth in vertical direction. The 
flow vector at cell is defined as 

^ i /-.,./E*(^;/;X)). 



Heuristic definition of local temperature for a one 
dimensional system 

In order to generalize Makse and Kurchan's definitions to mea- 
sure local temperature, we propose a heuristic procedure. First, 
we subdivide the whole system into smaller cells. Each cell has 
an index i. Next, we define the transition probability of a parti- 
cle from the ith cell to the jth cell P(j, i). Suppose the spacing 
between cells is £, a random walker must travel the distance of 
£ in order to go from cell i to cell i + Thus the averaged time 
to until this arises is obtained as 



to 



£ 2 /2D. 



The escape probability of a particle from cell i is 

P out (i) = ^ P(i + fe, i; / = 0) oc 1/to = 2D/£ 2 . 

k=±l 

On the other hand, when a drift force / exists, 

to = t/fli, 
then the flow J(t) along the / direction is 

J(i) =£^2 fcP (* + k,i;f^0)(x £/to = f/Jt. 
fc=±i 



where P(i', i; j) is the transition probability from to (z', j) 
and P{i;j',j) is the transition probability from to (i, j') 
(See Fig. [TJi. 



P(0;1,0). 



P(-1,0;0) 
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P(1,0;0) 



Fig. 1. Schematic figure to explain the relationship among the tran- 
sition probability P, the flow vector J and the unit vector e± per- 
pendicular to the flow vector. In this example, the center cell is 
numbered as (0,0). Thin solid arrows indicate the transition prob- 
ability, bold solid arrow indicates the flow vector, and broken ar- 
row indicates the unit vector e± perpendicular to the flow vector, 
along which the force is applied later. The flow vector J of the 
center cell is £{P{1, 0; 0) - P(-l, 0; 0), P(0; 1, 0) - P(0; -1, 0)). 
The escape probability of a particle from the center cell p out is 
(P(1,0;0)+P(-1,0;0),P(0;1,0) + P(0:-1,0)). 
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Next we can define the escape probability of a particle from 
the cell as 



P° ut (i>r,f = Q)= E P(i + k,i;j;f = 0). 
Similarly, 



fe=±i 



P° ut (i,3-J = Q)= £ P(i;j + k,j;f = 0). 

k=±l 

In order to deal with vibrated beds, we have to take into 
consideration the case when J is not zero. We have to find the 
direction along which the flow is zero when / is equal to zero. 
Such a direction is perpendicular to the flow vector. Thus the 
unit vector along this direction is 



e .(i a. f-Q)- J_ (- J y(hrJ = 0) 



The escape probability of a particle from the ith cell along this 
direction is defined as 

PT\hi;f = 0) = I e ± (i,j;f = 0) ■ P out ^rJ = 0)1, 



where 



P ou *(i,j;/ = o) = 



Next we apply a small force fe± at cell This time, 

the flow along this axis is defined as 

J±{i,j;f? 0) = e ± (i,j; / = 0) • J(i, j; / ^ 0), 

which takes non-zero values when / ^ 0. Using these above, 
we get 

VJ 2 Jx(*,i;/^0) ' 

as a local temperature at (i, j). Tfj is obtained after taking the 
limit of / — * as before. 




Fig. 2. The flow vectors J(i,j). The flow is upward at the center re- 
gion and downward at the wall region. The solid rectangle (4 < i < 
8, 3 < j < 9) shows the region where the effective temperature will 
be computed. The left-bottom corner corresponds to = (1, 1) 

and increases right(up)wards to imax = 11 and jmax = 13. 



gravity force is applied. This means that in a naive sense, it is 
difficult to expect a homogeneous state over the whole system. 
Thus it will be interesting to see whether thermal equilibrium 
is kept or not using the above defined effective temperature. 

In orderto compute T$ .first we compute P(i±l,i;j; f — 
0) and P(t; j ± 1, j; f = 0) numerically. Then applying /, we 
compute P(i±l, i; j; f ^ 0) and P(i; j±l,j; f ^ 0). / should 
be as small as possible, but too small / cannot generate large 
enough J±(i,j; f ^ 0) to be estimated. Also if applying / to 
all the cells at once, global flow structures can be destroyed. 
Furthermore, we need measurements for several forces before 
taking / — > limit. Considering these requirements, we apply 
f /mg = 0.1, 0.2, 0.4, 0.5 to only a few cells at once (m is the 
particle mass). Then repeat the same calculation after chang- 
ing the set of cells to which / is applied. Figure shows the 
dependence of upon f /mg at cells with j = 3. Since they 
have clear linear dependence upon /, we extrapolate to / = 
using a linear least square fit in order to compute T^. Since 
the dependence of TV, upon / / mg is similar for other layers, 
we omit the remaining plots. 
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Numerical measurement of the effective granular 
temperature 

The DEM method is the same as in the previous study 1 7 1. The 
number of particles is 100, the diameter of a particle is 1, and 
the horizontal size of the vessel is 1 1 . The coefficient of resti- 
tution is e = 0.79 and t c = 0.04, with k = 3000 (spring 
constant) and r\ = 6 (dissipation constant). The strength of in- 
teractions with the walls and the bottom is two times larger than 
the above k and r\. The acceleration amplitude P = 1.51g (g 
is the gravitational acceleration, here it is taken to be 98, am- 
plitude is 4.11 and angular frequency is 6.20). Figure |2] shows 
the J(i,j) averaged over sequential 4 x 10 8 snapshots, where 
£ = 1. This means N x = = 1, .., 11). N y can be taken as 
large as possible, but in the present study j = 1,.. .,N y (= 10). 
The (1,1) cell corresponds to the down-left corner, i and j in- 
creases rightwards and upwards, respectively. 

As can be seen easily, there is a convective motion over the 
whole system. The flow pattern is not uniform at all. Also the 
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Fig. 3. The dependence of Ty upon f/mg at j — 3. Solid lines in- 
dicate a linear extrapolation to / = using a least square fit. The 
numbers in parentheses indicate the targeted cells as (i, j). 



Figure |4] shows the spatial distribution of Zy . We have 
omitted the regions close to the wall and the bottom (4 < i < 8, 
3 < j < 9). Although If?* fluctuates from cell to cell, the 
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values seem to be close to each other (at least for i = 5, 6, 7). 
Thus surprisingly this suggests the system is more or less in the 
thermal equilibrium state. 
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eff 

Fig. 4. Spatial distribution of TR . Each row corresponds to a hori- 
zontal layer of cells. The whole region corresponds to the rectangle 
in Fig. [2] Missing values (i, j) — (4, 3), (7, 9) are omitted due to the 
large error. 



In order to confirm this conclusion, we propose the correc- 
tion equation of TfJ^. Why are these effective temperatures not 
constant? Does it mean the breakdown of thermal equilibrium? 
As mentioned above, should be measured along the axis 
without flow when / = 0. For this purpose, we selected the 
direction perpendicular to the flow direction J. However when 
measuring the effective temperature, a particle can be drifted to 
the neighboring cells where the flow along drift direction is not 
zero. In other words, in order that our formulation is satisfied 
well, the flow in each cell has to be as parallel to each other as 
it is in Makse and Kurchan's study. Of course, this requirement 
is not satisfied in the vibrated bed of powder at all. 

In order to remove the effect that comes from the fact that 
the flow vectors are not parallel to each other, we propose here 
the heuristic equation, 



1 + V2 



k 



sin 



-3) 



where To is the unperturbed true constant effective tempera- 
ture, and the summation is taken over the eight neighboring 
cells, k. Ok is the angle between the flow vector and the di- 
rection of / at cell k (Fig. The supposed meaning of this 
equation is as follows. sin(f5 fc — ?) is the measure of the flow 
along / direction relative to neighboring cells. It is propor- 
tional to the difference between the neighboring flow vectors 
if we ignore the difference of absolute values of the vectors as 
higher order corrections. Thus if it is negative(positive), par- 
ticles flow in the same(opposite) direction of /. This appar- 
ently increases(decreases) the drift caused by /, i.e. the de- 
nominator of the definition of the effective temperature, there- 
fore the measured temperature will be smaller (larger) than the 
true temperature To. In order to express this tendency we add 
sm(9k — §) to the right hand side. 1 + \[2 is just a numeric 



parameter that accounts for the square lattice. The factor 1 rep- 
resents distance between neighboring cells and v2 does that 
between next neighboring cells. It does not mean anything spe- 
cial at the present study. Figure|6]shows 



To = 



T eff (l 



V2) 



l + £ fc sin(0 fc -§) 



It takes an almost constant value for all cells. Considering that 
there are no fitting parameters, our heuristic argument turns out 
to be very good. Thus we can conclude that the vibrated bed of 
powder is surely in a thermal equilibrium state. 
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Fig. 5. Schematic figure to show the meaning of the correction equa- 

eff 

tion of the effective temperature Tfj . Solid arrows indicate the flow 
vectors at each cell. When we compute the effective temperature at the 
center cell, we apply a small force / (dotted arrow) to the center cell. 
However, if the flow vectors at surrounding cells are not parallel to that 
of the center cell, these accelerate or decelerate the drift caused by /. 
This occurs when the angle 6t between the flow vector of surrounding 
cell and the force direction differs from n/2. 
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Fig. 6. The corrected temperature as a function of the effective tem- 
eff 

perature . In spite of the large fluctuation of the effective temper- 
ature, the corrected temperature takes an almost constant value. 



Does segregation break the thermal equilibrium? 

In the previous section, the effective temperature calculated 
from the heuristic argument indicates that the vibrated bed of 
powder is in the thermal equilibrium state. Although the ki- 
netic property, i.e. flow, is non-uniform in the vibrated bed of 
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powder, particles are the same all over the vessel. Thus it is in- 
teresting to see if different particles can also be in the thermal 
equilibrium state or not. Actually, Makse and Kurchan 1 5 1 have 
shown that particles with distinct sizes have the same effective 
temperature. However, their system is kinetically uniform. 

In this section, we deal with the vibrated bed of powder 
with different particles. There are many ways to introduce dif- 
ferences among particles. The size difference is the most fa- 
mous in the vibrated bed of powder, because it causes segrega- 
tion. However, in our cell approach, different sizes may cause 
some difficulty because we have employed cells as large as 
each particle. Different particle size requires different size of 
cells or another approach. 

Let us remind readers that the mixture of glass beads and 
lead beads causes segregation 18 191 . This phenomenon has al- 
ready been reproduced by our model 1 10 1 1 1 when we use two 
kinds of particles with a distinct coefficient of restitution. Thus 
the segregated state by the introduction of the difference of the 
coefficient of restitution will provide us the opportunity to see 
that the thermal equilibrium is satisfied even if the system is 
spatially non-uniform in both kinetics and material properties. 
The coefficient of restitution used is e = 0.67 (t c = 0.04) [k — 
3200,7? = 10] and 0.88 (t c = 0.067) [k = 1200, 77 = 2]. 
The strength of interactions with walls and bottom is two times 
larger. The acceleration amplitude r is taken to be 1.62g. 

From the numerical point of view, a difference from the uni- 
form system occurs only at collisions between different parti- 
cles. When this occurs, k and 77 is taken to be the averaged value 
of these two. Everything else is the same as the previous nu- 
merical simulations. Figure0shows a snapshot and convective 
flow patterns. As can be seen easily, two kinds of particles are 
segregated and convection occurs in each layer independently. 

For the calculation of convective flow, we do not distin- 
guish different particles. The effective temperature is measured 
at each cell as above. Figure[8]shows the dependence of temper- 
ature at each cell upon /. Again linearity with / is very good 
and we get extrapolated values using the least square fits. Fig 
|9]shows the spatial distribution of temperature. Clearly spatial 
distribution is much more scattered than that in Fig. |4] This 
will be because of inhomogeneity of the system. However, ap- 
plying corrections to this as done in the previous section, we 
get almost constant values of the effective temperature (Fig. 
1 1 Oi . although they are a little more scattered than in the uni- 
form case. Thus it suggests that the thermal equilibrium state is 
not destroyed even if the system is not uniform. Even when the 
drastic phenomena like segregation occur, the thermal equilib- 
rium state can still be maintained. 



Discussion 

Although we have found a thermal equilibrium state using the 
temperature defined by us, there are still some uncertain points 
in this study. First, the heuristic arguments presented are not 
substantiated by any first-principles argumentation and the suc- 
cess of the study is judged based on the good uniformity of 
the effective temperature. Actually our finding that we have 
obtained a thermal equilibrium state does not always justify 
our heuristic argument of the temperature. However, even if 
we cannot justify our results completely, we believe that our 
findings cannot be accidental. It was very difficulty to find a 




Fig. 7. (a) A snapshot of a segregated bed. Open circles represent 
particles with the smaller coefficient of restitution. The two kinds of 
particles are segregated, (b) Flow vectors of the segregated bed. Con- 
vective motions are localized at each segregated layer. The solid rect- 
angle shows the region (3 < i < 9, 2 < j < 9) where the effective 
temperature will be computed. 
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Fig. 8. The dependence of Tij upon f /mg at j — 2. Solid lines in- 
dicate linear extrapolation to / = using the least square fit. The 
numbers in parentheses indicate the targeted cells as (i, j). 



phenomenological temperature definition which characterizes 
the thermal equilibrium state in the dynamical state of granular 
matter. Furthermore, our definition is based on the foregoing 
studies 1 5 6|, thus it has some physical bases although it is not 
very complete. Our findings can be a start point of understand- 
ing general dynamical properties of granular matter. 

Second, the numerical model studied here appears to be 
slightly artificial (two-dimensional, relatively few particles, cells 
of the same size as the particles, making it impossible to study 
systems composed of non-monodisperse particles) and there- 
fore the results obtained with this model are not as convincing 
as they might have been if a different, less restricted model 
had been employed. The reason why we had to employ such 
a model is mainly because of the computational limitation. In 
spite of the apparent easiness of our simulations, it is extremely 
time consuming. We need very long time simulations in order 
to suppress the fluctuations. This is because we need the ab- 
solute values of fluctuations, not mean values. Compared with 
getting mean values, measuring fluctuations quantitatively is 
much more time consuming. Furthermore, we have to repeat 
numerical simulations because we need extrapolations and the 
number of points which we can compute at the same time is 
only a few, as noted above. Of course, it is much more suitable 
to reproduce these results by more realistic numerical simula- 
tions, but it is beyond the scope of this paper. Validating them 
using more suitable systems will be another future issue. 
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Fig. 9. Spatial distribution of . Each row corresponds to each hor- 
izontal layer of cells. The whole region corresponds to the rectangular 
in Fig. |Zl Missing values ( = (3, 5), (4, 4), (8, 4), (9, 5)) are 

omitted due to the large error. 
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The effective temperature 

Fig. 10. The corrected temperature as a function of the effective tem- 
eff 

perature T% . In spite of the large fluctuation of the effective tem- 
perature, the corrected temperature takes almost constant value. The 
horizontal line indicates the mean value of the corrected temperature. 
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Summary and Conclusion 

In this study of a granular system, we have introduced heuris- 
tic procedures to measure the effective temperature (corrected) 
proposed by Makse and Kurchan 1 5 1 . In the vibrated bed, ther- 
mal equilibrium is maintained even if spatial uniformity is ki- 
netically destroyed by convective motion. Also it is still main- 
tained even when segregation occurs in the system. Since the 
effective temperature is introduced only phenomenologically, it 
is important to understand this phenomenon by statistical me- 
chanics. However, since clearly energy is not conserved, how 
to construct statistical mechanics in granular matter is unclear. 
Although we do not have any clear pictures at the moment, re- 
cent findings |12| strongly suggest that statistical mechanics 
is apparently valid even for the single particle system without 
collisions. Probably, there is some robustness in the apparent 
thermal equilibrium system of the dynamics of granular mat- 
ter. 



